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Abstract 

This paper deals with the theory of structure, stability, robustness, and stabilization for an 
appealing class of nonlinear systems which arises in the analysis of chemical networks. The results 
given here extend, but are also heavily based upon, certain previous work by Feinberg, Horn, and 
Jackson, of which a self-contained and streamlined exposition is included. The theoretical conclusions 
are illustrated through an application to the kinetic proofreading model proposed by McKeithan for 
T-cell receptor signal transduction. 
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1 Introduction 

This work was originally motivated by the study of the foUowing system of first-order ordinary differential 
equations: 

Co = fci (T* - Eto^O (M* - E^=oC^) - (fc-1,0 + fcp,o) Co 

Ci = fcp,i— lCi-1 — (fc— l,i + fcp,i) Ci 
Cm — kp,N~lC M — 1 l.AfC'jv 

where the subscripted /c's, as well as M* and T*, are arbitrary positive constants; the C^'s are nonnegative 
functions of time t, and dots indicate derivatives with respect to t. These equations arise in immunology, 
and describe a possible mechanism, due to McKeithan, cf. , that may explain the selectivity of T- 
cell interactions (some more details are provided below). McKeithan analyzed the equilibria of these 
equations, which represent steady-state regimes, mostly under the simplifying assumptions that kp^i = kp 
and k-i^i = fc_i for some fixed kp and During Carla Wofsy's series of talks on the topic, a 

number of questions arose: what can be said about the structure of equilibria? are the equilibria stable? 
Stability is a natural mathematical requirement, with clear biological significance. 

Under the simplifying assumptions mentioned above, it is not difficult to see that there is only one 
equihbrium, whose coordinates depend, of course, on the given constants. We will see that in fact even 

*Eduardo Sontag is with the Department of Mathematics, Rutgers University, New Brunswick, NJ 08903, USA. E-mail: 
sontagahilbert.rutgers.edu, URL; http://www.niath.rutgers.edu/ sontag This work was supported in part by US Air 
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in the general case of the constants being distinct, there is a unique equilibrium. Far more interestingly, 
however, we will prove that the equilibrium is globally asymptotically stable: every solution, for any 
choice of initial (nonnegative) values Ci(0)'s, converges to the unique equilibrium. This conclusion rules 
out, in particular, periodic orbits and, of course, chaotic behaviors, and shows the "determinism" of the 
process described by McKeithan. Moreover, we will also establish the robustness of stability with respect 
to a quantifiable class of perturbations in the dynamics. 

Going further, and approaching now the equations from the point of view of a control theorist, one 
may pose questions of a very different nature as well, questions of design. The equilibria depend, in 
particular, on the value of the constants T* and M*, which represent total concentrations (of intermedi- 
ate complexes plus T-cell receptors, and intermediates plus peptide-major histocompatibility complex, 
respectively). If these concentrations are allowed to vary, and are seen as additional state variables, 
multiple attracting equilibria exist. So one may ask: if we are allowed to manipulate some of the vari- 
ables, can we change the equilibria at will (and preserve stability)? How many variables need to be 
manipulated? Such questions might eventually impact approaches to therapy and rational drug design. 
We will provide a global feedback stabilization result as an answer. 

As it turns out, the system of interest can be viewed as a very special instance of a very large class 
of nonlinear systems, for which the above-mentioned results can be established in general. In this paper, 
we describe a theory of structure, stability, and stabilization for that general class of systems. Thus, this 
paper can be read totally independently of the above-mentioned motivating example from immunology. 
By providing a general theory, one expects that other applications of the theorems given here will be 
possible. 

The class of systems which we consider is, basically, that of "deficiency zero chemical reaction net- 
works with mass-action kinetics (and one linkage class)" in the language of the beautiful and powerful 
theory developed by Feinberg, Horn, and Jackson, cf. ^, |o[ |ll|, |2[ As a matter of fact, our stability 
(but not the robustness nor the feedback) results are basically contained in that previous work. The 
existence and uniqueness of equilibria, and local asymptotic stability, are already proved in the papers 
by Feinberg, Horn, and Jackson, and we claim no originality whatsoever in that regard. And although 
the global stability results may not be readily apparent from a casual reading of that literature, our 
proofs of them consist basically of a careful repackaging of the discussion found in Feinberg's paper j|] . 
(One should point out that an alternative approach to global stability, which would apply to a somewhat 
smaller class of systems, could use ideas from which, in turn, was based on |2^. Also somewhat 
related, but far more restrictive, are the results described in the book |l4[ , which apply to systems which 
satisiy a "conservation of energy" constraint.) 

The reliance upon the Feinberg/Horn/Jackson theory notwithstanding, we provide here a totally 
self-contained and streamlined exposition even of those results that are known for these systems (deep 
theories for other classes of systems have also been developed by these authors). We do this in order 
to present things in a terminology and formalism more standard in control and dynamical systems, and 
also because it is not obvious how to put together the necessary pieces from the various sources in a way 
that allows us to refer to them efficiently. In any case, it is our hope that this exposition will serve to 
make a wider audience in the dynamical systems and control theory communities aware of their work. 

In addition to the stability results, we also prove a robustness and a feedback stabilization result, 
and we show the existence of a global change of coordinates, which brings the systems in question into 
a canonical form which exhibits a particularly elegant structure. We develop a formalism, and present 
explicit estimates for stability margins, with a view towards further theoretical developments. Indeed, 
in a recent follow-up joint paper ||6[ with Madalena Chaves, we have been able to derive, using these 
techniques, various input-to-state stability results and, in particular, a design of globally convergent 
Luenberger-like observers for systems of the type studied here. As new results, we also include partial 



generalizations to non-mass action kinetics (cf. Remark 5.6 and the stability section) 



There is an extensive literature regarding applications of control theory to chemical engineering, 
and, in particular, to bioreactors. Among many references, one may mention the textbook Q and the 
more recent survey articles Q and and the references given there, as well as the well-known theory 
of compartmental systems (see e.g. pSj). More specifically, results based on Feinberg-Horn- Jackson 
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theory have appeared in the control hterature, see M. Other reactor biocontrol work deals with 
reachability and controllability issues, see for instance A preliminary version of this paper was 
posted electronically in p^ . 

The organization of this paper is as follows. In the rest of this Introduction, we motivate the formalism 
and results to follow by reviewing some basic facts concerning chemical reaction networks and working 
through an example. This section can be skipped with no loss of continuity by readers interested in the 
mathematical developments. The paper starts in Section |^, where we introduce the class of dynamical 
systems being studied, and the main theorems are also stated. In Section |[ we specialize to McKeithan's 
system, interpreting the results in that special case. Section || discusses some basic coordinatization facts, 
as well as some useful alternative system descriptions. Section H deals with the main proofs regarding 
interior equilibria. Section |6| with boundary equilibria. Section ^with technical invariance results, and, 
finally. Section ^ with the proofs of the stability theorems. 
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1.1 Chemical Networks 

Let us motivate the formalism and results to follow by reviewing some basic facts concerning chemical 
reaction networks. We will restrict attention to "mass action kinetics" ; however, it is important to remark 
that variations such as Michaelis-Menten reactions, obtained through singular perturbation analysis when 
starting from mass-action models, are routinely used to model many enzymatic reactions in biology. 

In chemical network studies, one analyzes systems of differential equations which describe the time- 
evolution of the concentrations xi{t), . . . , Xn{t) of n given "chemical species" Pi, . . . , P„; the Pi's might 
denote anything from small molecules to large complexes. The equations are derived from a consideration 
of the reactions that are known to occur among the substances Pi , perhaps helped by other molecules 
which are not explicitly considered in the equations (such as catalysts or energy sources). 

As an extremely simple illustration, suppose that each molecule of a certain species Pi can react with 
one molecule of P2 to produce a molecule of P3, and that, conversely, each molecule of P3 may dissociate 
(through a different process, typically with different time constants), into Pi and P2. This is indicated 
graphically by 

Pi + P2 ^ P3 and P3 ^ Pi + P2 , 

or just by 

Pi + P2 ^ P3 • 

An example of this reaction is provided by the synthesis of ethyl tert-butyl ether (P3 = CeHi^O) from 
isobutene (Pi = C^Hg,) and ethanol (P2 = C^H^O). For studies of the ethyl tert-butyl ether reaction in 
the control theory literature, see and [^. Assuming that the reactor is well-mixed, particle-collision 
theories or quantum-mechanical potential energy methods are often used to justify the statement that 
the probability of such a reaction occurring, in a small time interval around time t, is proportional to 
the product Xi{t)x2{t) of the concentrations at time t, and to the length of the interval, that is to say, 
to the probability of two molecules in this bimolecular reaction colliding by virtue of being "in the same 
place at the same time" . Since we gain a single molecule of P3 for each such reaction, we arrive at a 
formula for the rate of increase of the concentration of P3 due to the first reaction Pi + P2 ^ P3: 

±3 = fcia;ia;2 . (1) 

In addition, as one molecule of Pi and P2 each are eliminated at the same rate, we also have the following 
two differential equations: 

Xi ~ —kiXiX2 , X2 = —kiXiX2 ■ (2) 
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Here, ki is a suitable constant of proportionality, the "reaction rate constant" , which is often taken to 
be (Arrhenius law) proportional to the Boltzmann factor and T is the temperature. One also 

writes graphically: 

Pi + P2 ^ ■ 

We also assumed for our example that there is a dissociation reaction, that is, 

P3 -Pi + P2 

where k2 is another rate constant, and this reaction is also modeled by rate equations: P3 decays at a rate 
proportional to its concentration, ±3 = —k2X3, and each of xi and X2 grow at this rate. Incorporating 
these into the previous equations gives the final set of differential equations: 

ii = -kiXiX2 + ^2X3 

±2 = -kiXiX2 + k2Xz (3) 
±3 = kiXiX2 - k2X3, 

which describe the evolution of all the concentrations Xi{t). 

A very convenient and systematic formalism to describe the complete system of equations is as 
follows. The entire reaction is represented by a graph, whose nodes are the "complexes" which appear 
in the reactions, such as Pi + P2 and P3 in the example given above, and whose edges are labeled by 
the reaction rate constants. So, in the example shown above, there is an edge labeled fci (where fci is 
an actual positive number) starting at the node in the graph corresponding to Pi + P2 , and pointing 
to the node corresponding to P3, and there is likewise an edge labeled k2 from the second node to the 
first. We associate to this graph its 2 x 2 incidence (connectivity) matrix matrix A = {fly }, listing all 
the edge labels (for instance, 021 = to indicate a reaction with rate constant fci, from the first node 
to the second node). More generally, the size of A is m x m, if there are a total of m complexes. One 
of the hypotheses to be made is that the matrix A is irreducible, meaning that the graph is strongly 
connected, that is, there is some path, typically through several stages of intermediate reactions, linking 
any two given complexes (in fact, results hold under a weaker, block-irreducibility, property). Next, one 
introduces a set of column n- vectors 61, . . . , 6m, one for each complex (in our simple example, m — 2). 
This is done by specifying the contributions from each type of molecule. For example. Pi + P2 gives rise 
to the vector 61 = (1, 1, 0)', and P3 to the vector 62 = (0, 0, 1)'. 

The "mass action" dynamics are then summarized by the system 

m m 

* = EE a,,x\^^x\^^ ...x\-^{b,~b,) (4) 

i=i 3=1 

(where each bi £ M" has entries bu, . . . , bni), which the reader may easily verify reduces to (||) in the 
above example, for which 

-a oO.-(;).-(;). 

Note that, for instance, the term corresponding to i = 2 and j = 1 in the summation in (^ gives us 

(-kiXiX2 \ 
-fci2;ia;2 
kiXiX2 ) 

which provides precisely the contributions to ii, ±2, and ±3 represented by (|l|) and (||). (Higher-order 
polynomial equations may result as well. For example, suppose that each molecule of species Pi can 
react with four molecules of P2 to produce two molecules of P3, that is, P1+AP2 2P3. This would give 
us equations Xi = —kxix\, ±2 — —4:kxiX2, and 2:3 = 2fca;ia;|, and the vectors bi and 62 would become 
&i = (1,4,0)', 62^(0,0,2)'.) 
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A fundamental role is played by the linear subspace D of M" which is spanned by all the differences 
bi — bj. This is the stoichiometric subspace associated to the reaction, and each intersection between 
a parallel translate of T> and the positive orthant is called a class, or more properly a stoichiometric 
compatibility class ( "stoicheion" — (Gk.) element). The significance of classes is that, since i e D, 
trajectories remain in classes, that is, classes are positive-time invariant manifolds for the dynamical 
system. (The positive orthant is itself forward invariant, as is easily shown.) In the example discussed 
above, D is the line spanned by 61—62 = (1, 1, —1)'. This line can also be described as the set of solutions 
of xi + X3 = X2 + X3 = 0, so each class is given by the positive points in {xi + x^ = ci,X2 + X3 = C2}, 
for different constants ci and 02- Of course, it is clear from the equations (||) that d{xi + xz)/dt = 
d{x2 + X3)/dt = along solutions. 

One of the basic facts about the systems studied here is that there is a unique equilibrium in each 
class, and, under mild conditions, this equilibrium is globally asymptotically stable with respect to 
positive initial conditions. Continuing with the above example, the set of possible equilibria consists of 
the points in the hyperbolic paraboloid k2X3 — kiXiX2 = 0. Let's now take for simplicity ki ^ k2 = I 
and analyze the stability of the equilibrium a;" = (1, 1, 1)'. As the ray (1, 1, 1)' + A(l, 1, —1) is forward 
invariant, where A ranges over the interval (—1, 1) (the intersection of the corresponding line and the 
positive orthant), we may parameterize motions by A. One obtains the scalar differential equation 
A = — 3A — A^, which has, as claimed, the point A = as an asymptotically stable state with a basin of 
attraction which includes all of (—1, 1). 

If we start at a point which is not in the above line, the equilibrium approached will not be a;°. If 
this equilibrium is desirable, we may want to design a feedback law to drive the solutions from every 
other (positive orthant) initial state into cc". Suppose that we can control the inflows and outflows of, 
let us say. Pi and P3. This situation is represented by the following control system: 

Xi = -XiX2 + X3+ Ui, X2 = ~XiX2 + X3, ±3^X1X2 - 2:3 + "2- 

Stabilization of x^ ~ (1,1,1)' can be achieved for instance by a simple feedback linearization, taking 
Ml :— X1X2 — X3 — a;i + 1 and U2 —x\X2 + 1. Then, x\{t) 1 and x^it) — > 1 as t — > cx), and thus, by 
a cascade argument, we see that also X2(t') — > 1 as t 00, as wanted. We will provide a general result 
on stabilization (but using linear feedback), as well as results on robustness of stability and on global 
decompositions of dynamics, for systems (^. These results will be illustrated again for the above simple 
example, now viewed as a particular case of McKeithan's system, in Section ^. 

2 Definitions and Statements of Main Results 

Some standard notations to be used are: 

• R>o (resp., M+) = nonnegative (resp., positive) real numbers 

• M" (resp., M™^™) — n-column vectors (resp., m x m matrices) with entries on M+; similarly for 
M>o 

• Rq = boundary of R^O' ^'^^ '-'^ vectors x £ ]R>q such that xt = for at least one i £ {1, . . . , n} 

• x' — transpose of vector or matrix x 

• \x\ = Euclidean norm of vector in R" 

• (x, z) = x' z, inner product of two vectors 

• D-L = {x\{x,z) = OVz G D}. 
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Although we develop the theory for a somewhat wider class of systems, we wish to emphasize that all 
the results to be given are valid, in particular, for the following general class of systems evolving on K"q: 

rn rn 

X = aij ■ ■ ■ Xj^' {bi — bj) 

i=i j=i 

(and, as appropriate, for those systems obtained by adding control inputs, as discussed below). Each 
column vector b£ 6 R" has entries bu, . . . , bni, which are nonnegative integers, and the a^j's are non- 
negative numbers. The systems defined in this fashion are described by polynomial dynamics. The only 
assumptions required in order for the results to hold are that the bi^s be linearly independent and that 
the m X m matrix A — (oij) must be irreducible. Recall that this means that (/ + A)™~^ G R™^™ 
or, equivalently, that the incidence graph G{A) is strongly connected (where G{A) is the graph whose 
nodes are the integers {1, . . . ,m} and for which there is an edge j — > i, i 7^ j, if and only if a.y > 0). 
This assumption amounts to a "weak reversibility" property in the application to chemical reactors, as 
discussed briefly in Section |[ and is crucial to the validity of the results. 

We now describe the underlying dynamics of the systems to be studied, and leave for later the 
introduction of additional terms in order to model the possibility of control actions. Our systems are 
parametrized by two matrices A and B with nonnegative entries, as well as a collection of nonnegative 
functions 9i, i = 1, . . . , n, and have the following general form: 

m m 

i = f{^) = EE ^li^lf^'H^^)"'' ■ ■ ■ OniXn)'"' {h - b,) (5) 
i=l j=l 

where b( denotes the £-th column of B (notice that the diagonal entries of A are irrelevant, since 
bi — bi =0). Several restrictions on A, B, and the OiS are imposed below. The powers are interpreted 
as follows, for any r, c > 0: r" = 1, 0^ = if c > 0, and = e^'"'' if r > and c > 0. 

The main motivating example, arising from mass-action kinetics in chemistry, is obtained when 
9i{y) = \y\ for all i and i? is a matrix whose entries are nonnegative integers (so, for nonnegative vectors 
X, we have polynomial equations). This is the case mentioned in the first paragraph, and will be referred 
to as the "standard setup" in this paper. 

The hypotheses on the 6*2 's. A, and B are as follows. Each map 

Oi : [0, 00) 

is locally Lipschitz, has Oi{0) — 0, satisfies \\n9i{y)\ dy < 00, and its restriction to E>o is strictly 
increasing and onto. We suppose that 

A = (ay ) e M^o^" is irreducible (6) 

and, for B, that: 

each entry of B is either or > 1 (7) 

(this last hypothesis insures that fix) in (|) is a locally Lipschitz vector field, so we have uniqueness of 
solutions for the differential equation), 

B = (5i, . . . , 5,„) G R^^"" has rank m (8) 

(so, its columns bi are linearly independent), and 

no row of B vanishes. (9) 

This last hypothesis is made mainly for convenience. Observe that if some row, let us say the fc-th one, 
were zero, then the same dynamics would be obtained if all entries in row k are replaced by "1" (since 
the differences b^i — b^j are still zero) and one restricts the dynamics to those states satisfying Xk = c, 
where c is the positive number such that 9k{c) — 1. 
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From now on, we assume that all systems considered satisfy the above assumptions. 

Our study will focus on those solutions of (||) which evolve in the nonnegative orthant IR>q. Recall 
that a subset S C M" is said to be forward invariant with respect to the differential equation x = f{x) 
provided that each solution x(-) with a;(0) 6 S has the property that x(t) £ S for all positive t in the 
domain of definition of x(-). It is routine to show (cf. Section |^) that the nonnegative and positive 
orthants are forward invariant: 

Lemma 2.1 Both R"q and R" are forward-invariant sets with respect to the system (||). 

(These properties are simple consequences of the fact that, because of the assumptions made, the 
fc-th co mpon ent of a solution of (||) will satisfy Xk{t) > whenever Xk{t) = 0.) We will also show in 
Section ^.2.4 that there are no finite explosion times: 



Lemma 2.2 For each ^ G M"g there is a (unique) solution x{-) of (||) with a;(0) = ^, defined for all 
t>0. 

In order to state concisely the main results for systems (^, we need to introduce a few additional 
objects. The subspace 

D :— span {bi — bj, i ^ j} — span {bi — b2, . ■ . ,bi — 6™} (10) 

can be seen as a distribution in the tangent space of R"; it has dimension m — 1 because adding bi to 
the last-shown generating set gives the column space of the rank-m matrix B. For each vector p G M", 
we may also consider the parallel translate of D that passes through p, i.e. p + T) = {p + d, d^ 2)}. A 
set S which arises as an intersection of such an affine subspace with the nonnegative orthant: 

S = (p + D)f|M|o 

(for some p, without loss of generality in ]R>o) will be referred to as a class. If S intersects the positive 
orthant R!j:, we say that S* is a positive class. The significance of classes is given by the fact that any 
solution x(-) of (||) must satisfy 

xit)-xiO) = / f{x{s))ds = y y n{t){b.-b,) eD, (11) 



pf mm 

/ f{x{s))ds = yj2^{t)ib,-b,) eD, 

•^0 ^=l , = 1 



K{t) = Jg fly 6*1 (xi («))''« 6*2 (a;2(s))''2j . . . 6'„(x„(s))''"j ds, so x{t) £ a;(0) + T) for all t. In particular: 

Lemma 2.3 Each class is forward invariant for (||). □ 

The introduction of control action, through additional feedback loops, may be used in order to 
overcome the constraints imposed by (0). In order to formulate our control-theoretic results, we will 
suppose that r external inputs ue, can be used to independently infiuence each of € state coordinates. 
In other words, we will also consider the following control system associated to the basic open- loop 
model 

r 

X = f{x) + ^ uiCki (12) 
i=i 

(/ is as in (H)), where r is a positive integer, ei, . . . , e„ are the n canonical basis vectors in R", and 
fci, . . . , fcr are r distinct elements of {1, . . . , n}. 

Of course, inputs might influence the system in manners other than through independent action 
on some coordinates, which represent what are sometimes (especially in compartniental models) called 
inflow- controlled systems, see [y. Bilinear control action, for instance, is arguably more interesting in 
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reaction systems. We will only study the above class of control systems, but research is ongoing on 
generalizations to other formulations. 

We will first state the main results for autonomous systems (H), and later we state a stabilization 
result for (|l|). 

We denote by E (respectively, i?+ or Eq) the set of nonnegative (respectively, positive or boundary) 
equilibria of (|5|), i.e. the set of states x G M"g (respectively, e M" or G Kq), such that f(x) = 0. Of 
course, E is the disjoint union of E^ and Eq. 

Theorem 1 Consider any system under the stated assumptions. For every maximal solution of 
with x{0) G ]R"q, it holds that x{t) E as t —> +00. 



This will be proved in Section 8.2.2. The invariance of classes (which are contained in subspaces of 
dimension m — 1 < n) precludes asymptotic stability of equilibria of (|5|). The appropriate concept is 
that of asymptotic stability relative to a class. We say that an equilibrium x G S* is asymptotically stable 
relative to a class S if it is (a) stable relative to S (for each e > 0, there is some 6 > such that, for all 
solutions x(-), |a;(0) — x\ < 6 and x(0) G 5* imply \x{t) — x\ < e for all t > 0) and (b) locally attractive 
relative to S (for some e > 0, if |a;(0) — x\ < s and a;(0) G S then x{t) ^ x as t ^ +00). We say 
that a; G 5* is globally asymptotically stable relative to a class S if it is stable relative to S and globally 
attractive relative to S {x{t) — > x for all solutions with a;(0) G S). The main results are as follows; the 



first part is shown in Section p.l| , and the remaining two in Section |8|. 

Theorem 2 Consider any system under the stated assumptions. Fix any positive class S. 

a. There is a unique equilibrium xs G 

b. The equilibrium xg is asymptotically stable relative to S . 

c. The equilibrium xs is globally asymptotically stable relative to S if and only if Sf]EQ — 9. 



Example 2.4 The following two trivial examples may help in understanding the above theorems. In 
both cases we take n = m = 2 and 9i(i/) = \y\ for i = 1,2. The first example has 



The system (||) is (for nonnegative states): 

±1 = (xi - l)xiX2 
X2 = 

and thus Eq = Rq = {x \ X1X2 = 0} and E+ = {x\xi = 1,X2 > 0}. The positive classes are the sets 
S = Sr = {x \ xi > 0, X2 = r}, for each r > 0, and for each such S ^ Sr, xs = (1, r)' is asymptotically 
stable with domain of attraction {x\xi > 0, X2 = r}. See Figure |l|. Each class Sr has a second 
equilibrium (0, r)', but this second equilibrium is in the boundary, so there is no contradiction with part 
a of Theorem |^. Regarding Theorem 0, observe that every trajectory either converges to an interior 
equilibrium (1, r)' or it is itself a trajectory consisting of an equilibrium (and hence also converges to E, 
in a trivial sense). 

The second example has 
We now have (for nonnegative states) the following linear system: 

±1 = X2- Xi 
X2 = Xi- X2 
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1 r 
Figure 1: Equilibria for example 2.4, dark lines indicate equilibria 



and thus Eq — {0} and = {a; | a;i ^ X2 > 0}. The positive classes contain no boundary equilibria; 
they are the sets S = Sr = {x M>o \xi — X2 = r}, for each r > (see Figure |l|). For each such S = Sr, 
xs = {r,ry is globally asymptotically stable. □ 



Feedback stabilization 

The invariance of classes precludes the existence of global (not merely relative to a class) attractors for 
the uncontrolled system Now suppose that we wish to find a feedback law which forces all the 
solutions of the closed-loop system obtained from the controlled system (|l2|) under this feedback to 
converge, as t — » +oo, to a specific equilibrium x G M!ji. It turns out that the obvious negative feedback 
solution, namely to use inputs proportional to the errors Xk — Xk, achieves the goal of global stabilization 
to X, provided that enough inputs are used (r > n — m + I, which is obviously necessary, since solutions 
of (^2|) evolve in linear subspaces of dimension m — 1 + r) and that the fc^ 's are appropriately chosen. 
We next state a result in that regard. For each j G {1, . . . , m} we consider the set 



8, := {k I bkj > 0} (13) 



which is nonempty, by (| 



Theorem 3 Let r = n — m + 1, and suppose that 71, ... ,7^ are arbitrary positive real numbers, and 
ki, . . . ,kr G {1, . . . , rt} are such that 

D + span {efe,,...,efcj = (14) 

and 

8j C {fci,...,fcj (15) 
for some j G {1, . . . , m}. Pick any equilibrium x G E-\-. Then, all maximal solutions of 

r 

X = f{x) + ^ 7^ {xke - Xki.) ekt (16) 

with a;(0) G R"^ are defined for all t > and remain in M>o; cind x is a globally asymptotically stable 
equilibrium of (Uqj. 

The invariance statement is proved in Section |^. The global stability statement is proved in Sec- 
tion ^.2.1| . 



To study the closed- loop system (|l^) , and also to be able to formulate a result concerning robustness 
of the stability properties described by Theorems |l| and |[ we will study "positive perturbations" of the 
basic uncontrolled system model (||). These are described by equations as follows: 

i = f*{x) = f{x)+9{x) (17) 
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where / is as in (||) (and all the stated assumptions hold) , and 5 is a locally Lipschitz vector field on R" 
for which the following property holds: 

(VxeM|o) (Vfce {!,..., n}) [xfc = ^ gfc(x) > 0] (18) 

where gkix) denotes the fc-th coordinate of g{x). Property ( |l^ ) is the most natural assumption which 
guarantees the forward-invariance of ]R"q (since it will imply that ik{t) > whenever Xk{t) = 0). Of 
course, any results established for arbitrary systems (|l^) will be also true for systems (^) (take g = 0). 

The feedback system (^ is of this type. To see that (^ holds, note that either gk{x) — ^i{xkf —Xk^) 
(if k = some kg) or gk{x) — 0, and, in the first case, Xk ^ Q implies gk{x) — ^yiXk^ > 0. Moreover, it 
satisfies the following strengthening of (|l8|): 

(3je{l,...,m}) (VxeRlo) (VfceSj) [xk ^ ^ gk{x) > Q] (19) 

(note the strict inequality, in contrast to ([T^)). Indeed, pick j as in (|l5|), k §j, and x so that x^ — 0. 
Then, for each k € §j, property ( p^ ) guarantees that gk{x) = ^iXk^ > 0. 



Robustness 

A different specialization of the general form ( pT| ) allows the study of robustness with respect to pertur- 
bations which preserve classes. The corresponding systems are obtained by adding vector fields which 
lie pointwise in the span of the hi — 6j's. We suppose given a collection of locally Lipschitz functions 
Ay : M" M>o (i, j £ {1, . . . , m}) such that 

(Vi,je{l,...,m}) (Va;GM|o) (VfceSj) [xfc = ^ Ay(a;) = 0] (20) 

and, using these, define the system: 

m m 

i = /(2^) + EE^^^(^) (^*-^^) • (21) 

Observe that property ( ^ ) implies that this system is of the general type ([l^). Indeed, the only possible 
negative signs for gk{x) = J2 ^iji^)(fiki — bkj) can arise from the terms of the form Aij{x)bkj with 
bkj 7^ (i.e. k e Sj), but these vanish, because of (^), when Xk = 0. Of course, systems ^ are a 
subclass of ( ^Tl ) (take all A^j = 0). The main robustness result will be as follows. 

Theorem 4 For each positive class S there exists a continuous function Sg : S'P|M" R>oi with 
6s{x) > if and only if x ^ xs, such that, for any collection {Aij} such that 

rn rn 

EE^^^-(^)' ^ (22) 

i=i j=i 



for all x G 5P|M" , the following properties hold for the system ^2]\): 

1. BothWyQ andW^, and the class S, are forward-invariant. 

2. For each ^ G S* there is a (unique) solution x{-) with x{0) — ^, defined for all t > 0. 

3. The equilibrium xs is asymptotically stable relative to S . 

4- The equilibrium xs is globally asymptotically stable relative to S if and only if Sf]EQ = 9■ 
ln other words, the claim is that the conclusions previously stated for systems (^ are preserved by 
perturbations, as long as the magnitude constraint ( p^ ) is satisfied. The result is not in itself surprising; 
the main interest is in its proof, which will be largely constructive, offering an explicit formula for the 
function ^5 and, moreover, the fact that 63 will depend nicely on the class S (in a sense that will be 
clear). The invariance of S is clear from the form (pl|), and the (also easy) invariance of orthants is 
proved in Section 0. Parts 2-4 are proved in Section 3.2.3. 
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Regularity 

For all results to be stated next, we suppose a fixed system (||) has been given. 

We say that a function defined on an open subset of K" is of class C'^, k = 1,2, .. . ,00, w if it is 
/c-times continuously differentiable (for k = 1, . . . , 00) or real-analytic (for k = lu). We will refer to the 
following hypothesis, for each such k: 

(Hfe) each 9i restricted to R+ is of class C'^, and 9^{y) > for all y > 0. 

Obviously, this hypothesis is satisfied with k — uj when 9i is the identity map for y > 0, as is the case in 
the standard setup leading to polynomial systems. 



We will show, in Section 5.2 



Theorem 5 If hypothesis (H^) holds, then £'+ is an embedded submanifold o/R" , -diffeomorphic to 

j^n — m+l 

Observe that, in the standard setup, / is a polynomial vector field, so £'+ = {x G R" | f{x) = 0} and 
hence is an algebraic subset. 

Theorem H states that there is a unique interior equilibrium x in each positive class. This equilibrium 
depends nicely (smoothly, analytically) on the class, provided that the 9iS be regular enough. In order 
to make this statement precise, we introduce the map tt : K" — » R" which assigns, to each x G R" , the 
unique interior equilibrium xs in the class S which contains x. The next result is proved in Section ^.l| . 

Theorem 6 If hypothesis (Hk) holds, then tt is of class C'^. 



Note that, once that we know that an embedded submanifold, there is no ambiguity in the above 
statement: if we view tt as a map tt : M" it is also of class C''. 

The following open subset of R" : 

n := M.I + D (23) 

is the union of all those parallel translates of D which intersect the positive orthant. It includes all 
positive classes. Assuming again hypothesis (Hfc), it turns out that one may always find a change of 
variables which transforms TZ onto R", and in particular E+ onto the set {(Xi, X2) | = 0} and positive 
classes into subsets of sets of the form {(Xi, X2) | X2 — c} for constants c, and transforms the dynamics 
on TZ into the following form: 

Xi ~ Fi{Xi,X2) 
X2 ^ 

where Xi and X2 represent blocks of variables according to the decomposition R" — R'"^^ x M"^™+^. 
Since points of arc equilibria, this implies that i^i(0, X2) = 0. To be precise, we prove in Section ^.3| : 

Theorem 7 Assume that hypothesis (Hk) holds. Then, there exists a diffeomorphism $:??.—> 
R''"-i X R^-^+i such that, denoting 

F{X) = ci>,($-i(X))/($-i(^)) 
and writing F = {Fi,F2) and $ ~ ($1, $2) in block form, 

1. X — z Cz D if and only if $2(2;) — ^2iz), 

2. $1(2;) = if and only if x G E^, and hence F{Q,X2) — for all X2, and 
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3. F2{X) = for allX. 



The preceding results provide steps in proving this one, but conversely, Theorem ^ implies Theorem ^, 
since E+ is transformed into a coordinate space. In fact, it also establishes that E+ transversally 
intersects each class S. Also Theorem ^ implies Theorem ^, since the map that selects the element in 

in the same class as a given point amounts to a coordinate projection, let us call it tt, onto the last 
n — m + 1 variables under the change of variables given by Moreover, let O be the set of no nneg ative 
points that are not boundary equilibria of (||), that is, O = R"q \ Eq. We will prove (Corollary^ that 
each X G O belongs to some positive class; thus, the mapping tt extends to O. Since TZ includes every 
positive class, it includes O, and the extension of tt to O transforms into the restriction of vr to ^{O). 
Thus, we have the following consequence of Theorem |7[ 

Corollary 2.5 If hypothesis (H^) holds, then the extension of tt to O is also of class . □ 



3 Mass- Action Kinetics, and McKeithan's System 

As mentioned in the Introduction, the results explained in Theorems |l| and ^ are basically theorems 
for what are called mass-action networks of zero deficiency, and are given (implicitly in the case of 
global stabihty) in [p| l ioi |ll], |l§ . (There is one global stability result stated explicitly in the above 
papers, namely in Wt\. The statement would be that every trajectory which starts in the positive 
orthant must converge to the interior equilibrium in the corresponding class. However, this question is 
still open, since the suggested proof used the implication that if a positive definite function V of states 
has a negative derivative along trajectories while away from an equilibrium, then global stability of the 
equilibrium follows. This implication is, in general, false, without the assumption of some sort of radial 
unboundcdness, a property which definitely does not hold in the context in which it is being applied; 
see [|l^ for a retraction of that proof.) 

The ass ump tions of irreducibility of A and full rank of B (both of which may be relaxed somewhat, 
cf. Remark |5.5| ) are key ones. They serve to rule out periodic (or even chaotic) behaviors which may 
otherwise arise in chemical networks such as the Belousov-Zhabotinsky reaction or Prigogine-Lefever's 
"Brusselator" (for which see e.g. 0). 



3.1 Kinetic Proofreading in T-Cell Signaling 

The equations with which we started represent the dynamics of the "kinetic proofreading" model pro- 
posed by McKeithan in in order to describe how a chain of modifications of the T-cell receptor 
complex, via tyrosine phosphorylation and other reactions, may give rise to both increased sensitivity 
and selectivity of response. Let us introduce two additional variables T{t) and M{t), which represent the 
concentrations of T-cell receptor (TCR) and a peptide-major histocompatibility complex (MHC). The 
constant ki is the association rate constant for the reaction which produces an initial ligand-receptor 
complex Co from TCR's and MHC's. The quantities Ci{t) represent concentrations of various intermedi- 
ate complexes, and McKeithan postulates that recognition signals are determined by the concentrations 
of the final complex Cw The constants kp^i are the rate constants for each of the steps of phosphorylation 
or other intermediate modifications, and the constants are dissociation rates. 

Global stability of a unique equilibrium will be deduced from Part c in Theorem |^ when we view the 
equations as those of an appropriate system (^ restricted to a suitable class (which is determined by 
the constants M* and T*). The complete reaction network is represented graphically in Figure || (this 
is the same as Figure 1 in JT^ , except that we do not make the simplifying assumption of equal rates) . 
In other words, we have a system of form (||), where we use x = (T, M, Cq, . . . , Cjv)' as a state, and take 
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Figure 2: McKeithan's network 
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flij = 0. The corresponding set of differential equations is: 



t = -kiTM + ^k_i^,C, 



M 
Co 



i=0 



i=0 



kiTM - (fc_i,o + fcp,o) Co 



C^ 



Notice that D = {x | T + Co + . . . 
form S' = 5"'^, intersections with 



Cjv — ^p,JV — iCjv — 1 ^— l,JvC*jv ■ 

. . + C„ = ^/ + Co + . . . + C„ — 0}, and the positive classes are of the 
"q of the afhne planes 

T + Co + . . . + C„ = a , M + Co + . . . + C„ /3 , 



with a > and /3 > 0. The original system is nothing else than the restriction of the dynamics to 
the class determined by a = T* and /3 = M* . Thus, the conclusions will follow from Theorem ^ as 
soon as we prove that S* f] £Jo = for any positive class S. Pick any a; G Kg and any positive class 
8°"'^ . According to Proposition |6.3| , it will be enough to find some j £ {1, . . . , m} with the property 
that Xk 7^ for all fc £ Sj. Here, Si = {1, 2} and %j = {j + 1} for j = 2, . . . , m. If the property is not 
satisfied for some j G {2, . . . , m}, then Ci ~ for all i. But in this case, the equations for 5'"''^ give that 
T — a > Q and also M = /3 > 0, so = 1 can be used. In conclusion, x ^ Eq, and hence Part c of the 
theorem applies. 

Let us compute equilibria explicitly. Setting right-hand sides to zero gives CjY — (/bp^^_i//c_i jY)CjY_i, 
and recursively using Ci — + kp^i))Ci-i we may express all C^'s as multiples of Cq. We 

also have Co = [fci/(fc_i,o + kp^Q)]TM. Thus 



{{a, j3, K,Qaj3, Kia/3, . . . , Kf^a[3) | a, /? > 0} 



where the k^'s are rational functions of the constants defining the system (these are the equilibria studied 
in [p^). It is obvious in this example that we obtained a two-dimensional (n — to + 1 = 2) nonsingular 
algebraic subvariety of M" . 
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We next discuss the application of Theorem ^ to this example. As r = 2, we consider feedback laws 
of the type jiixk^ — Xk-^) + "f2ixk2 ~ 2:^2): where 7^ > and fci, k2 satisfy (14) and (p^). Since V is 
spanned by ei + 62 — 63, . . . , ei + 62 — e„, this means that dij) is satisfied provided that 1 e {ki, or 
2 g {ki,k2}, and this is sufficient to guarantee that also (|15D satisfied. In other words, stabilization to 
any desired equilibrium is possible provided that some pair of reactants, including at least one of T or 
M, can be manipulated by a controller. 

To illustrate some of our constructions in a comparatively trivial case, let us now specialize even 
further, taking just = and all constants equal to one. (This is mathematically the same as the 



example treated in Section |1.1| .) Writing x, y, z instead of T, M, Co, we have the system 

X ~ —xy + z , y = —xy + z , z = xy — z 

and E-f- = {(x, y, z) \ z — xy = 0} is a hyperbolic paraboloid (intersected with the main orthant). 

This is obviously a nonsingular algebraic set, but it is instructive to see it a lso as the image of the 
diffeomorphism M : constructed in the proof of Theorem || in Section ^.2[ Using x = (1,1, 1)', 

we have the formula M(ri, r2, ra) = (e''^ , e^^, e'^^)' for (ri, r2, r^) G D^. Here, D is the span of (1, 1, — 1)', 
i.e. {(ri,r2,r3) | ri + = r2 + = 0}, and = 7'2, ^3) ks = ?'i + ?'2}, so it is clear that 
M{D^) = E+. 

Note, incidentally, that there are boundary equilibria as well: Eq = {(a;, 0, 0) | a; > 0} 2/1 0) | y > 

0}, but none of these is in a positive class (since both x + z and y + z are positive constants on positive 
classes). 

The most natural (in the motivating application, anyway) initial state is one in which z(0) = 0. We 



now compute the value of the map tt considered in Corollary 2.5 at such a point (0, 2/0,2:0). We must 
find the positive-orthant intersection of {z = xy} with the line Lxo,yo = {x + z = xq, y + z = y^}. 
Equivalently, we need to solve {z — xo){z — j/o) — z — subject to the constraint < z < minjxo, ?/o}j 
which guarantees that x > and y > as well as z > 0. There is exactly one such solution, and it is 
the smallest of the two solutions, since the graphs of fi{z) — z and /2(z) — {z ~ xq){z — j/o) intersect 
at precisely one point in the interval (0, minjxo, yo}), so we take the negative sign in the quadratic 
formula: 

1 



7r3(xo, 2/0,0) = - xo + I/O + 1 - V (2^0 - yo^ + 2(xo + ya) + 1 
(and corresponding values for tti and 112), which is, indeed, a real-analytic (hypothesis Hj^) function. 

Finally, let us find a decomposition as in Theorem]^. We let Xi :— z and X2 — (u, v) = {x + z,y + z). 
Thus X2 = along all solutions, and Xi — z — Fi{z,u,v) = {u — z){v — z) ~ z. This does not yet 
have Fi{Q,u^v) — for all so one needs a further change of variables to bring the equilibrium to 
zero, which can be achieved by translating z ^-^ z := z — 7r3(0, u,w). Note that the set of points 
with z = xy, gets mapped into z = 0, and the positive orthant gets mapped into the set given by the 
inequalities m > 0, u > 0, and z G (— 7r3(0, u, w), min{u, w} — 713(0, u, w)). 

4 Some Preliminary Facts 

The equations (|^) (and their perturbed form (|l7|)) have a considerable amount of structure, and various 
useful properties are reflected in alternative expressions for the system equations. 

4.1 Other Expressions for the System Equations 

For each i = 1, . . . , n, let us introduce the maps 

Pi{y) ■■= In (2/) 
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(with pi{0) = — oo); note that liniyx^o Piiv) = ^co and \pi{y)\ dy < oo. Furthermore, the restriction 
of Pi to is locally Lipschitz, strictly increasing, and onto R. For any positive integer n, we let 

p : K" ^ [-00, cx))" : x i-> (pi(a;i), . . . ,p„(a;„))' 

(we do not write "/9„" to emphasize the dependence on n, because n will be clear from the context). 
Then (|) can also be written as 

m m 

i = /(^) = EE"''^-^^'^'"^'^^ - • (24) 

The expression "e^''^'^^^^^" is interpreted in accordance with the conventions made for powers: if a; is a 
vector and fc e {1, . . . , n} is an index such that x^ — and b^j > 0, then e^''^^'''^^''^ = 0, consistently 
with e~°° — 0, and thus also 



e 



{bj,p{x)) _ ^bij pi{xi) ^b2j P2{x2) ^ ^ ^^Kj P,zix„) _ Q 



but, if bkj = 0, then we have e'''=3'''=(^'') = 1. 

Another useful way of rewriting (|^) is as follows. We write fk for the fc-th coordinate of / (i.e., the 
coordinates Xk of solutions x satisfy Xk = fk{x)). The terms in the sums defining fk can be collected 
into two disjoint sets: those that do not involve a product containing Ok{xk), for which bkj = 0, and 
those which do involve 9k{xk). The latter, by assumption (^, have bkj > 1, so we can factor 9k{xk) 
from ei{xi)''^ie2[x2f''' ■ . .6l„(a;„)^"^ and there remains a locally Lipschitz product. In other words, we 
can introduce, for each k G {1, . . . , n}, these two locally Lipschitz functions: 

Mx) -.^ E (E«y ^1(2^1)'^^ ^2(0:2)''^^... efe(a;fe)'"=^-^..0„(a;„)^- (25) 

and 

Pkix) := (E ^i(2;i)''^^2(a;2)''^^ . . .0„(x„)''- , (26) 

where Jk.\ '■— {j | bkj > 1} and J^^o ■= {j \ bkj = 0}. (Note that, in terms of the sets introduced in (13), 
k G §j if and only if j e Jk,i-) Note that (3k{x) > for all x. More generally, for perturbed systems (17), 
we let 

/3fe := Pk+ gk- 

In terms of these functions, 

fk{x) = ak{x)6k{xk) + Pk{x) , 
fkix) = akix)0kixk)+f3kix) 

(where fk=fk+ gk)- If a; S lR>o and k are such that Xk — 0, then property ( p8[ ) says that gk{x) > 0. 
In particular, since 0fc(O) = 0, 

Xk = Q ^ fk{x) = I3k{x) > and /,!(x) = pkix) > (27) 
so the vector fields / and /* = f + g always point towards the nonnegative orthant, on the boundary 

4.2 A Coordinatization Property 

When we apply the following Lemma, we will always take D = T), but we can state the result in more 
generality. Actually, we may get an even more general result by dropping the assumption that the maps 
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6i are onto. We will only assume, for the next result, that each i = 1,. . .,n, Oi : M ^ [0,(7^) (where 
< (Ti < oo) is locally Lipschitz, has 9i{0) — 0, satisfies |ln6'j;(?/)l dy < oo, and its restriction to R>o is 
strictly increasing and onto [0, di), but we are not now requiring ct^ = +oo. The reason for this relaxation 
is to allow consideration of functions like 9{x) — which arise in Michaelis-Menten kinetics. We will 
impose, instead, another condition, which we describe next. Let, as earlier, Piijj) := InOiijj), and let 
yOj := IncTi (infinity if (Xj = oo). Notice that , seen as the inverse of the restriction of pi to is a 
strictly increasing map from (— oo,/?;) onto M+. Thus, for any given constant p, p~^{s) > p + 1 for all 
s > to, some to G (— oo, pi), which implies that, for L{t) — P^^{s) ds—pt, its derivative {dL/dt){t) > 1 
for all t G (to, Pi). Under the general assumption pi — oo, we have that, for any p: 



lim 



I p:r\s)ds-pt = +00 (28) 



for any finite a < pi. Instead of assuming pi = oo, we will merely ask that (28) should hold. (Example: 



x) = fefj gives p ^(s) = so L{t) = -fcln(l - e*) - pt ^ +oo as t ^ = p.) 

Lemma 4.1 Let D be any subspace of R". For each p, q in , there exists a unique x = (p{p, q) € R" 
such that: 

x-p e D (29) 

and 

p{x) - p{q) e . (30) 
Furthermore, if hypothesis (Hfc) holds, then the map p,q i-^ v{Pi l) is also of class C'^. 



Proof. We fix p, q as in the statement, and start by introducing the following mapping, for each i G 
{l,...,n}: 

Li{t) := / p~^[s)ds - Pit 

Jq 

defined for t ^ Ci := {—oo, pi — pi{qi)). By assumption (|2^), Li{t) increases to infinity ast y pi — pi{qi). 
Also, Li(t) — > +00 as t ^ —oo, because p^^ is nonnegative and pi > 0. Thus, Li is proper, that is, 
{t \Li{t) < v} is compact for each v. 

Now we take the (continuously differentiable) function 

n 

Q{y) := 5]L,(2/,) 

i=l 

thought of as a function oi y E C := Y[7=i ^i- Observe that is in the open set C, since pi — Pi{qi) > 
of all i by definition of pi. This function is also proper, because 

n 

{yeC\ Q{y) <w} C e A I L,{t) < w - [n - !)£} 

1=1 

where £ is any common lower bound for the functions Li. Restricted to D^f^C, Q is still proper, 
so it attains a minimum at some point y G p| C, which is a nonempty (since belongs to it) 
relatively open subset of D^. In particular, y must be a critical point of Q restricted to D^f]£, so 
(VQ(y))' G (D^)-^ = D, which means that 

Pl^{yi+ Pl{qi))-Pl,---,Pn^{yn+ Pn{qn))-Pn) ^D- (31) 

Finally, pick x G such that p{x) = y + p{q). Such an x exists because we can solve the equations 
Pi{xi) = yi + Pi{qi) for each i, because yi < pi — Pi{qi) since y E C Then p{x) — p{q) G D-^ by definition, 
and (|3l| ) gives also x — p E D. 
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Finally, we show uniqueness. Suppose that there is a second z G M" so that z — p G Z) and 
p(z) — p{q) G D-^ . This implies that x — z D and p{x) — p{z) G D^. Since each pi is an increasing 
function, we have that, for any two distinct numbers a, b, (a — b){pi{a) — Pi{b)) > 0. So 

n 

'^{xi - Zi){pi{xi) - pi{zi)) = {x - z,p{x) - p{z)) = 

i=l 

implies x — z. 

It follows that (p{p,q) = a; is a well-defined mapping from M" x M" into M!f.. Suppose now that 
hypothesis (Rk) holds, and therefore also each pi has positive derivative and is C*"' for positive arguments. 
Let be an n X (m — 1) matrix whose columns are a basis of D (for instance, when D = we may 
take the columns hi — b2, ■ ■ ■ ,bi — bm), and let ^ be an n x (n — m + 1) matrix whose columns form a 
basis of D^. Consider the map F : MJ} R" given by 



F{x) 



V'x 
W'p{x) 



Observe that F is of class C'^. Denote by J{x) the Jacobian of F evaluated at an a; G K." . 

We claim that J{x) is nonsingular, for any x. The transpose of J{x) can be expressed as a block 
matrix: J^x)' — {V,TW), where T ~ diag{p[(xi), . . . , p'„{xn)) is a symmetric positive definite matrix. 
As the columns of TW are linearly independent, as are the columns of V, it will suffice to prove that 
the column spans of V and of TW intersect only at zero. Since the column spaces of V and W are 
orthogonal, this fact follows from the following observation: if T is self-adjoint and positive definite, then 
TM^Pl W^-'- = {0} for any subspace W of R". (Proof: factor T = R'R, with R nonsingular, and suppose 
that Tw G TWC]W^. For any WQeW,0= {Tw,wo) = {Rw,Rwo). Thus, Rw G RW f]{RW)^ = {0}, 
so Rw = 0, which implies w = 0.) 

To conclude, we note that, given p and q, x ^ viPi <?) is the imique solution of 

Gix,p,q) = F{x) ~ (^^'P^^^ = 0. (32) 

Thus, the Implicit Function Theorem gives that is class C'^, since ^{x,p, q) = J{x) has rank n at all 
x,p, q and G is of class C'^. I 

Remark 4.2 Note that, since G{(p{p, q),P; q) = in (^2|), we know that 

Ji^ip,q))^ip,q) + (^') = (33) 



dp"-"^'^' \ 

for all p,q. Thus {d(p/dp){p,q) = —J{(p{p,q))^^{V,Oy, and we conclude that dip/dp has constant rank, 
equal to n — TO -|- 1. □ 



The following quantity measures deviations relative to D . Let us define, for each x, z G R" : 

[{bi - bj,pix) - p{z))^ 

i=l 3 = 1 



S{x, ^) := E E - - '°^(^))) • (34) 



Remark 4.3 Note that 6{x, z) = if and only if p{x) — p{z) G D-^, since D is spanned by the differences 
bi — bj. Also, note that, since x — x = OgT>,x satisfies both ( p9| ) and ( ^ ) (which uniquely characterize 
fix, z)) if and only if p{x) — p{z) G . To summarize: 

X = ip{x, z) ^ 5{x, z) = Q ^ p{x) - p{z) G T)^ (35) 

for any x, z G R" . □ 
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Remark 4.4 We could also have defined a smaller, but basically equivalent, sum using only the gener- 
ating differences bi — bi, i — 1, . . . ,m — 1, but the above definition for d seems more natural. Moreover, 
note that if we let 



m— 1 



S{x,z) := ^ {{bi - bi,p{x) - p{z))f + ^ {{ve,x- z)f 



where the Vi constitute a basis of , then the uniqueness part of Lemma 4.1, applied with = 2), 
gives that 5{x, z) = if and only \i x = z. (Because x = (p{x, x), and 5_{x, z) — implies z = (p{x, x).) □ 



Remark 4.5 Recall the definition (^) of TZ. We can extend <y9 to a map TZ x M" M" by defining 
^p{p, q) '.— (p{p + d, q) for any p E TZ and q G R" , where d G D is any element of T> with the property 
that p + d G . This definition is valid because, for any di and d2 in 2) such that p + di and p + d2 
are in R" , we have (p{p + di,q) = (p(p + d2,q), since x = f{p + di,q) is uniquely determined by the two 
conditions 6{x, q) — and x — {p + d^) G D, and the second condition is equivalent to a; — (p + 1^2) G D. 
Moreover, if hypothesis (H^) holds, then the extension of 1^9 to 7?. x R" , which we denote again by (p, is of 
class C'^. Indeed, given any po G TZ, and a do G 1* so that po + do G R" , it holds that p + do £ R" for all 
p in some neighborhood of po , so the extension of ip is obtained (on that neighborhood) by composition 



of the original p with the translation {p, q) {p + do, q), and the former is of class C by Lemma 4.1 



Remark 4.6 If hypothesis (Hi) holds, then, for some continuous d{-) > 0, and for all x, z € 



z eT> 



5{x,z) > d{z) \x — z\^ + o ^|a; — 



(36) 



Indeed, S{x,z) > J2iL2i{^i ~ bi,p{x) — p{z))Y , and p{x) — p{z) + Q{x — z) + o{\x — z\), where Q = 
Q{z) = diag(pi(zi), . . . ,pj^(z„)) is the Jacobian of p evaluated at z (nonsingular, by hypothesis (Hi)), 
so 5{x, z) > J2iL2ii^i ~ bi,x — z)qY + o{\x — z|^), where (ry, S)q := (ry, QS) is a nondegenerate inner 
product on 2). This gives (|3^), because the bi — biS constitute a basis of D. (Or, more concretely: 
^{x, z) > YZ2 WAt + o{\x - z\^), where ^ := R.{x - z), {R,{b2 - 61), . . . , R.{bm ~ 61)), and 

Rz is a symmetric positive definite matrix, continuous on z, with R^ — Q(z). So, we may use that 
\Wz^\ > do{z) 1^1, where do{z) is the smallest singular value of restricted to RzT), which is nonzero 



since W^^ = implies ^ G (-R^D)-^, together with \^\ > di{z)\x- 
eigenvalue of Rz, to obtain again an estimate (36).) 



where di{z) is the smallest 



□ 



5 Interior Equilibria 

It is convenient to also express the dynamics (^ in matrix terms. Letting 



A = A-diag^a,a,...,^( 



we write: 



where Qb is the mapping 



4=1 



/ - l^i^l «il ai2 
^21 - Z^^/2 



V 



O-ml 



0'm2 



DTI , ju>m 



X = fix) = BAesix) 



p(x)) Jb„,,p{x)) 

..... 



air. 

a2r. 



Em . 



(37) 



18 



obtained as the composition of the maps x i— > p{x)^ z ^ B'z, and y > (e^^, . . . ,e^™)'. In particular, 
since ranki? — m and each pi maps M+ onto M, 

the restriction 65 : M" ^ M" is onto. (38) 

Note that 

Qsix) = [Oiixif'^e^ix^f-' . . . eniXn)'"',- ■ ■ , ^1 (o^i )"-" 02 (2:2 )'^" ■ • ■ (x„ )"'"") ' ■ 

For any two x,x E M" , it holds that: 

(3 k > 0) ei3(x) = Kei3(x) ^ /9(x) -p(i) e T)^ (39) 

(recall the definition ( |l0[ ) of 2)). To see this, denote y := Qsix), y :— Qb{x)- If y = Ky, then, with 
k :— InK, In^j = fc + Inj/j, for each j = 1,. . . ,m. Thus {bj,p{x)) = k + {bj,p{x)), which implies that 
{bj,p{x) — p{x)) — k for all j, and therefore 

{hj -h,,p{x) - p{x)) ViJ. 

Conversely, if this holds, we may define k := (hi, p{x) — p{x)), k :— e'' , and reverse all implications. 

We also note, using once again that B has full column rank, that f{x) = is equivalent to AQb{x) = 
0, that is: 

Lemma 5.1 A state x is an equilibrium if and only if 0_b(x) £ ker A. □ 

We now consider the matrix A. The row vector 1 — (1, . . . , 1) has the property that lA = (0, . . . , 0), 
so, in particular, A is singular. The following is a routine consequence of the Perron-Frobenius (or finite 
dimensional Krein-Rutman) Theorem. 

Lemma 5.2 There exists y E f] ker A so that {W^^ \ {0}) f] ker 1 = {ny, k > 0}. 

Proof. If y E IR.>Q is any eigenvector of A, corresponding to an eigenvalue A, it follows that = lAy = 
IXy = Aq, where q :— ly is a positive number (because y, being an eigenvector, is nonzero), and therefore 
necessarily A = 0. In other words, a nonnegative eigenvector can only be associated to the zero eigenvalue. 
Pick now any 7 > large enough such that all entries of A := A+jI are nonnegative. Since the incidence 
graph G{A) coincides with G{A), it follows that A is also irreducible. By the Perron-Frobenius Theorem, 
the spectral radius ct of A is positive and it is an eigenvalue of A of algebraic multiplicity one, with an 
associated positive eigenvector y E M™. Moreover, every nonnegative eigenvector y E M>q associated to 
cr is a positive multiple of y. As adding 7/ moves eigenvalues by 7 while preserving eigenvectors (that 
is, {A + jl)y = (A + 7)2/ is the same as Ay — Ay), y is a positive eigenvector of the original matrix A. 
It is necessarily in the kernel of A, since we already remarked that any nonnegative eigenvector must be 
associated to zero. Finally, if y is any other nonnegative eigenvector of A, and in particular any element 
of (IR>Q \ {0}) Pi ker A, then it is also a nonnegative eigenvector of A, and thus it must be a positive 
multiple of y, completing the proof. I 



Corollary 5.3 The set of positive equilibria is nonempty. Moreover, pick any fixed x E E^. Then, 
for any positive vector x E K" , the following equivalence holds: 



X E E< 



S{x, x) = 0. 



(40) 



Proof. By Lemma 5.2, there is some y E K.™ in ker^. By (|38|), there is some x E M!fi such that 85(2;) — y. 



In view of Lemma 5.1, x is an equilibrium 



Now fix any x E E^ and any x E M", and let y := Qb{x), y := Qsix). Suppose x E £'+. By 
Lemma 5.1, y E ker^l. By Lemma ^.2| , every two positive eigenvectors of A are multiples of each other, 
so there is some k E K+ such that y ~ ny. By (p9[), p{x) — p{x) E T)^. Conversely, if this holds, then, 
again by (|39|), y = ny. hence y is also an eigenvector of A, so by Lemma 5.1 we conclude x E I 
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Remark 5.4 Even if the assumption that each 0i\Ry.g is onto R>o is dropped, the second part of Corol- 
lary 5.3 is still valid, as the proof of the previous lemmas did not use ontoness. We are not assured that 
equilibria exist, but if there is any equilibrium x then (40) holds. □ 



Remark 5.5 It is possible to generalize many of the results to the case when A is not irreducible 
but, instead, there exists a permutation matrix P with the property that P' AP is a block matrix with 
irreducible blocks. Let us sketch this next. Assuming already such a reordering has taken place, the 
family of systems considered is as follows. The dynamics are described by 

i = fi{x) + . . . + fUx) 

where L is a positive integer (the "number of linkage classes"), and each fs has the form in Equation (||), 
for some matrices A = As and B = Bg. We suppose that each As is an irreducible nonnegative matrix 
of size ms x m^, and each Bs has size n x m^. We let m — mi + . . . + m^. Further, we also assume that 
the matrix B which is composed from the blocks Bg, namely B = (Bi, . . . ,Bl), satisfies the properties 
(0)' ®)^ ® required of _B. We let Ds be the subspace of R" spanned by the differences 6^ — &j of columns 
of Bs, for each s = 1, . . . , L, and now define 2) as the sum of the spaces D^. (Note that the column 
spaces of the S^'s intersect only at zero, because we are assuming that the columns of B are linearly 
independent; thus 2) is also a direct sum of the D^'s. The dimension of each is rus — 1, and D has 
dimension m — L.) Since the differences bi — bj of columns of different i?s's do not ever appear in the 
vector fields defining the system, the same argument as before shows that the cosets p + D (with the 
new definition of 25) are invariant. Classes are now defined using this D, and are also invariant. We can 
express the dynamics in the form (|37|), where A is formed as before, starting from the m x m matrix A 
that is obtained by using Ai, . . . , Al as diagonal blocks. 

Since each Ag is irreducible, we can find positive eigenvectors for the entire matrix A, and the 
possible such eigenvectors are of the form (kij/J^, . . . , K^y^)', where each ys is a positive eigenvector 
of As and the k^'s are positive numbers. The function Ob is still onto, so we obtain the existence of 
positive equilibria. Note that 9b (x) decomposes into blocks Qg, . . . , mapping into R>g respectively, 
and ( p9[ ) generafizes to: &%{x) = K&g{x) for some k > iff p{x) — p{x) G D^. Thus, uniqueness in 
each class holds just as before, since D"*- is the intersection of the 25^, s = I, . . . ,L, and also x is an 
equilibrium iS all fi{x) — 0. 

As a last comment along these lines, we remark that the assumption that B has full rank can also 
be slightly relaxed, as follows. Suppose that the column spaces of the B^'s intersect only at the origin, 
and the column space of some Bs spans an affine space of dimension rUs — 1, i.e., Vs has dimension 
nis — 1, but Bs itself has rank nis — 1 (instead of uis). Then, we may add a state variable to the system, 
which satisfies a differential equation Xq = 0, and extend the system in such a manner that Bs has rank 
nis'. this is equivalent to adding a row constantly equal to 1 to the matrix Bs (and to the remaining 
matrices as well). The original system appears as the restriction of the new system to the invariant 
subset consisting of those states whose last coordinate is equal to one. □ 



5.1 Proof of Part a in Theorem |2| and Theorem |^ 



Fix an X € and any p e R^fi. Applied with q = x (and D — D), Lemma LI says that x — (p{p, x) is 
uniquely characterized by (a) x £ p + 2) and (b) d{x,x) — 0, the second of which is in turn equivalent 
by (|4^) to X £ E^. In conclusion, x — (p{p, x) is the unique point in R" which is an equilibrium and lies 
in the same class as p. This proves Part a in Theorem ^ and, moreover, shows that 

Tr{p) = (p{p,x) (41) 

is the assignment of p into the unique positive equilibrium in its class. Furthermore, by Lemma |4.l| , 
ip{-,x) is of class provided hypothesis (H^) holds, which proves Theorem^. I 
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Remark 5.6 Theorem ^ is valid even under the weaker assumption that condition (|2|) holds, not 
necessarily assuming that each 6i is onto M>o. More precisely, i/ there is at least one equilibrium, then 
the conclusions hold. The stability parts will be remarked upon later, but regarding Part a, it is only 



Remark 5.4) 



necessary to note that Lemma 4.1 was proved under the weaker hypothesis, and that ([40D holds (cf. 



Now let us introduce the following condition, expressed in terms of the transpose of the vector 
1 = (1, . . . , 1) and the image of the transpose of B: 

1' G B'{Rl) . (42) 

(Observe that, if the system is "homogeneous" in the sense that the exponents in each term add to a fixed 
constant, then ( ^2|) holds, since homogeneity amounts to the requirement that i?'l„ = 1, where we denote 
by in the vector of length n consisting of I's.) We claim that if (|4|) holds (instead of unboundedness), 
then there exists at least one interior equilibrium. For this we may assume, without loss of generality, 
that (Ti = 1 for all i (if this were not the case, we simply rescale all the functions 9i, adjusting the A 
matrix accordingly; thus we rewrite the system equations as those of a system of the same form, for 
which (Ti = 1). It follows that the image of Qb is exp(i?'(— M" )), where "exp" means taking exponentials 
of each coordinate. Now take any y S M™ P|ker^ and let v S R™ be so that y — exp(ti). As i3'(M") is 
an open set, and using (|4^), there are A > and q G K" such that 1' — {1/X)v = B'q, or equivalently 
— Al + V — B'{—Xq) which means that e^^y — Qb{x) for some x G K". As e^^y G ker A, x is an 
equilibrium, as wanted. □ 



Remark 5.7 Recall Remark 4.5. After fixing an x G we may extend tt to 7?. by letting Tr{p) := 
ip{p,x) for any p € TZ. As (p(j>, x) := ip{p + d, x) for some d so that p + d G R" , and (p{p + d,x) G i?+ 
by (^l|), we have that tt{p) G for all p dTZ. Note also that n{p) = p ii p € since in that case p is 
itself the unique positive equilibrium in its class. Thus, 

p = n{p) ^ peE+ (43) 

for any p E TZ. Furthermore, as (fij), x) — {p + d) = ip[p + d, x) — {p + d) G © by definition of ip, we know 
that (p{p, x) E p + T>, so 

p - n{p) G D (44) 

for all p G TZ. We also have that 

7r(pi) ^ n{p2) Pi - P2 e 2) (45) 

for all pi,p2 G TZ, by ([44[) and because 7r(p2) is the unique positive equilibrium 'm p2 + D ~ pi + 'D when 



Pi — P2 G D. Finally, note that tt is of class C if hypothesis (Hfc) holds, once more by Remark 4.5. □ 



5.2 Proof of Theorem |5| 

Suppose that hypothesis (H^:) holds. Then the restriction of each pi to R+ is of class C*^, and so is 
its inverse, which is a map R — > R"*". With some abuse of notation, we will denote in this section the 
restriction of each pi to R+ , and more generally the associated vector function p, restricted to R" , with 
the same symbols. Note that p is a C'^-diffeomorphism (class C*^ and has an inverse which is also of class 
C''). Now we fix any equilibrium x G and in terms of it, define: 

M:R"^R!^ : y ^ THv + Pi^)) , 

which is a C'^-diffeomorphism since both y i— > y+p{x) and p~^ are C'^-diffeomorphisms. As an illustration, 
in the standard setup we would have Pi{r) = Inr and p~^{s) — e'*, so M{y) — (Sie^i, . . . ,a;„e^'')'. 

We claim that M(2)^) — E^. Once that this is shown. Theorem]^ will be established, since dif- 
feomorphisms map submanifolds into submanifolds, and has dimension n — (m — 1). Note that 
X G M('D^) if and only if p{x) = y + p{x) for some y G CD^, that is, if and only if p(x) — p{x) G 2)^, so 
Corollary |5.3| gives the equivalence with x G Ej^ . I 
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5.3 Proof of Theorem 



We will construct a preliminary transformation $ of 7?. into D x E+. Since 2) is a linear space of 
dimension to — 1, and since by Theorem || we know that is a manifold diffeomorphic to M"^™+^ , the 
result will then easily follow. We assume that hypothesis (Hfe) holds, and fix any x € E+. 



Recalling Remark 5.7, we define, for a; G 7?.: 

$(a;) = ($i(a;), $2(2;)) := (x - 7r(a;), 7r(a;)) 
(see Figure H). The function $ is of class as a map into M" x M" , and, because of 





2i\ 




^T) + x 














X — Tr{x)\ 









Figure 3: Change of coordinates 



x€E+ <^ $i(x) = ^ = (0,2:) (46) 

for all X ^ TZ. Also, property ( ^4| ) gives that <i>i(a;) G D for all x, and we have that ^2ix) G 7?+ (since tt 
assigns positive equilibria) . So the image of $ is contained in 2) x , which is an embedded submanifold 
of M" X M" because D is a linear subspace and using Theorem |5[ Thus, $ defines also a C*^ mapping 
into the submanifold T> x and we view it as such. 

Now consider the mapping 

* : D X M!^ ^ M" : {a,b)^a + b 

and its image 7^ = \E'(D x R^) = + CD. If we show that the restriction ^ of ^' to 2) x E+ is the inverse 
of $, then we will know that $ is a diffeomorphism. It is clear from the formulas that ^'(<I>(a;)) = x. 
Pick now any (a, 6) G 2) x E^, and let x = 'i'{a,b) = a + b. Thus, x — a = b <E E+, and (p{b, x) = b 
by ( |3^ ) and (pO|). By the definition of (/?, and thus tt, on TZ, we have that Tr(x) — ip{x, x) = (p{x — a, x),so 
^(x) — {x — (fix — a, x), ip{x — a, x)), which equals (a; — b,b) = (a, 6). Thus also 'I>(\I'(a, b)) = (a, b). 

In conclusion, we have a diffeomorphism 7?. — > 2) x , which can be composed with a diffeomorphism 
-01 X ■02 : 23 X K™^^ x (where -01 is a linear map) to yield the final diffeomorphism 

$ = ($i,$2) : 7^ M™"^ X M"-"'+i. We let F be the vector field obtained under this change of 
coordinates (namely, F{X) = <i>*(<i>^^(X))/($"^(X))) and write F = (^1,^2) in block form. Note that. 
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by (|4^), $i(a;) = if and only if x £ E+. That is, block vectors (Xi, X2) with Xi = correspond under 
the coordinate change to That is, every vector of the form {0,X2) is an equilibrium. This implies 
that F{0, X2) = for all X2, and in particular the same is true for the first block Fi. Also, by ( |45| ) we 
know that a; — z G 2) if and only if $2(2^) = ^2(2), i.e., if and only if $2(2;) = ^2{z)- 

Finally, we show that F2 = 0. It suffices to prove that X2{t) = <i>2(a;(t)) is constant, for an arbitrary 
solution x{t) of (||) taking values in TZ. This follows from the fact that classes are invariant, and hence 
the function ip{x{t), x) is constant, which gives from the definition that ^2{x{t)), and thus also ^2{x{t)), 
is constant. I 



6 Boundary Equilibria 

Fix any x = {xi, . . . , Xn) G ^>q- We wish to study the implications of some coordinate Xk vanishing. 
Recall the definition ( p^ ) of the sets Sj. We will use repeatedly the following fact, for any j e {1, . . . , to} 
and fc G {1, . . . , n}: 

Xfc = and fc G 8j => 9i{xi)'''' 92{x2p' ■ ■ ■ 6'„(a;„)''"^ = (47) 
which is obvious, since 9k{xk)''''^ vanishes when a;^ = and bkj 7^ 0. In particular, 

Xk = ^ 9i{x^f''92{x2p' . ..9M'-^bk, = (48) 
since either bkj = or ( ^7|) applies. We state results for arbitrary perturbed systems (|l^). (Note 



that ( [4 91) below is merely the particular case of (50) that arises for systems with 5 = 0, but it will be 



useful to have the conclusion for fk stated separately.) 

Lemma 6.1 Take any i, j G {1, . . . , m} such that aij ^ 0. Then, 

(V^ G §,) (a;, > 0) (VfcGS,)(a;fc>Oor/fe(a;)>0) (49) 

and 

(V£ G §j) {xi > 0) (Vfc G %i) {xk > or /*(x) > 0) . (50) 

Proof. Pick any k G Si, and assume that Xk = 0. The assumption "£ G §j ^ xi > 0" is equivalent to 
j e Jli ^ xi>Q. So, since Xk = 0, necessarily j G JkS)- By (|^), fk{x) = I3k{x), where Pk is as in (|2^), 
and the index j being considered does appear in the sum defining [3k- Moreover, for each £ G Sj, ^ 
by hypothesis, so 9i[xi)^^^ 92{x2)^'^^ . . . 0„(a;„)''"J > 0. On the other hand, k e §i means that bki > 0, 
and also ^ 0. Thus, the term involving this particular i and j in the sums defining [3k{x) is positive 
(and the remaining terms are nonnegative). We conclude that I3k{x) > 0, establishing (Ha). In addition. 



since we are assuming that Xk — 0, property (27) also says that fki^) = Pk{x) + gk{x)] as we proved 
that I3k{x) > 0, and gk{x) > (by assumption (|lS|)), also ( ^ ) holds. I 

We will denote by E* (respectively, E'^ or Eq) the set of nonnegative (respectively, positive or 
boundary) equilibria of (ll7|), i.e. the set of states x G K>o (respectively, G M" or G Rq), such that 
f*{x) = 0. When g = 0, these are the same as the sets E, E+, and Eq, and, more generally, for 
systems (|l7|), we still use E, etc., for the equilibria of the system (||) which has the same /. 

Proposition 6.2 For any system (0), and for an arbitrary x G IR>o, consider the following properties: 

1. xeE^. 

2. (VjG {!,..., to}) (3 fcG§,)xfc = 0. 

3. (Vj G {1,...,to}) ^i(xi)''i^^2(x2)'^^ ...0„(x„)^- =0. 
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4. X e Eq. 



Then 1 =^ 2 O 3 ^ 4. 



Proof. [1 2] Pick any x € Eq. If the second property is false, then there is some index j such that 
a;^ > for all £ ^ §j. We claim that for every index i, Xk > for all k G §i. Since IJ^ S^- = {1, . . . , n} 
(recall hypothesis (Q)), this will mean that Xk > for all fc, so x could not have been a boundary point, 
a contradiction. Let J = {j \ xg > 0\/ £ & Sj} and let / = {1, . . . , m} \ J. We know that J ^ and must 
prove that / = 0. Suppose by contradiction that / 7^ 0. Pick some i e / and j £ J such that a^- ^ 
(irreducibihty of A), and take any k £ Si. We claim that Xfc > 0. Suppose that this is not the case, i.e. 



Xk = 0. From ( pO[) in Lemma 6J, we conclude that /^(x) > 0. This contradicts the fact that x is in E*. 
In conclusion, Xk > for all fc £ S^, contradicting the fact that is/. 

[2 <^ 3] The product in 3 can vanish if and only if some bkj > and Xk — for the same fc, but this 
is precisely the condition in 2. 

[3 ^ 4] Since all terms in the definition (||) of f{x) vanish, it follows that x £ E. On the other hand, 
X must be a boundary point, since otherwise no product as in 3 could vanish. I 

For systems (^, all the properties are equivalent, since Eq — Eq. More generally, this stronger result 
holds for all systems of form : 



Proposition 6.3 For any system (^Tj) and for an arbitrary a; E IR>oi the properties in Proposition 6.2 
are equivalent. 



Proof. We prove that 4 => 2 => 1. If 4 holds, we apply Proposition 6.2 to the system (|^) which has the 



same /, to conclude that 2 holds. Assume now that 2 (and hence also 4) holds. Since we already know 
that f{x) — 0, all that needs to be verified is that Ay (a;) — for all Pick any Choose k G §j 
so that Xk = 0. Property (pO|) then gives the conclusion. I 



7 Invariance 

We start with an easy but key observation. 

Lemma 7.1 Suppose that x : [Q,t*] M"q is any solution of (p^). Then the following implication 
holds for any k = 1, . . . , n: 

Xk{0)>0 => Xk{t*)>0. 

Proof. Suppose that k is so that Xk{0) > 0. Let F : — ^ M be the (locally Lipschitz) function which 
coincides for t G [0, t*] and y £ K with 

Fit,y) fki^iit), ■ ■ ■ ^^k~i{t),y,Xk+i{t), . . . ,Xn{t)) 

and has F{t,y) = F{0,y) for i < and F(t,y) = F(t*,y) for t > t*. Note that F{t,0) > for aU t, 
because ( p7| ) says that /^(a;) = /3k{x) > when Xk = 0. For t £ [0,t*], the scalar function y{t) :— Xk{t) 
satisfies the scalar differential equation 

m = F{t,y{t)). 

We must prove that y never vanishes. For this, we let G{t,p) F{t,p) — F{t,0) and introduce the 
auxiliary initial value problem 

z{t) - G{t,z{t)), z(0) = y(0). 
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Since G is locally Lipschitz (in both variables, but just on z locally uniformly on t would be enough) and 
is an equilibrium of i = G{t, z), z{t) > for all t in its domain of definition. Moreover, we have that 
z = G{t, z) < F(t, z) for all t. By a standard comparison theorem, see e.g. [|l5|. Corollary 1.6.3, we know 
that z(t) < y{t) for all t in the common domain of definition of z and y. Since y(t*) is well-defined, z(t) 
remains bounded, and thus is defined as well for t — t* . So, y{t*) > z{t*) > 0. I 

Corollary 7.2 The set M" is forward invariant for (|T 

Proof. Consider any solution x : [0,T] R" of (|l7|), and suppose that a;(0) G M". We must prove 
that x{T) G R". Since x{0) is in the interior of M>q, the only way that the conclusion could fail 
is if x{t) € M.Q for some t G (0,T]. We assume that this happens and derive a contradiction. Let 
t* := minjt G [0, T] \ x{t) G MJJ} > 0. By minimality, for aU i, x,{t) > for aU t G [0,t*), and in 
particular Xi{t) G K"q for all t G [0, t*], and also there is some index k such that Xk{t*) = 0. But this 



contradicts the conclusion of Lemma 7.1 



The closure of an invariant set is also invariant, so: 
Corollary 7.3 The set R"q is forward invariant for 



Note that Corollaries 7.2 and 7.2 prove Lemma 2.1 as well as the invariance parts of Theorems ^ 
and ||. 

Lemma 7.4 Consider any solution x : [0, T] M" of ( p^ for which a;(0) G IR>o- Suppose that there is 
some jo G {1, . . . , m} such that 

(V^ G §jj [xe{0) > or /;(a;(0)) > 0] . 
Then, x{t) G M.^ for all t G (0, T]. 



Proof. We know that x{t) G K>oi by Corollary 7^, and we need to prove that Xk{t) > for all k and all 
t G (0,T]. Let 

/ ;= {ie {l,...,m}\xkit) > Vfc G §^,yt G (o,r]} . 

Since IJj §i = {1, . . . , n}, it will be enough to show that / = {1, . . . , to}. 

We start by remarking that jo G /, because, for each k G Sjg, either 2:^(0) > 0, and then Lemma [tt] , 
applied on any subinterval [0,i*], says that Xk{t) > for all t, or Xk(0) — /fe(a;(0)) > and a;fc(0) — 



imply that Xk{t) > for all t small enough, so also (again by Lemma 7.1) for all t. 

Suppose by way of contradiction that H := {1, . . . , to} \ / ^ 0. Pick some i G / and h € H so that 
o-hi 7^ (irreducibility of A). We will show that, for any given Iq G (0, T], and for any given k G Sh, 
Xk{to) > 0, and this will contradict h E H. 

Since i G /, Xi{to/2) > for all £ E Then, we can apply Lemma to obtain that Xk{to/2) > 
or f^{x(to/2)) > 0. As before, if Xk{to/2) > then via Lemma ^t] we conclude the positivity of Xk(t) 
for all t > to/2, and in particular of Xk{to). And if instead fk{x{to/2)) > 0, then Xk{to/'^) > and 



Xfc(to/2) > imply Xk{t) > for all t > to/2 near to/2, and so by Lemma 7.1 once again Xk{to) > 



For the special case of systems (^, and more generally perturbed systems (^ij), we conclude that 
every trajectory starting on the boundary Rq which is not an equilibrium must immediately enter the 
positive orthant: 



Corollary 7.5 Consider any solution x : [0,T] R'^g of 1^) for which a;(0) ^ E*. Then, x{t) G 
for aU i G (0,r]. 
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Proof. By Proposition 6.3 



a;(0) ^ Eq implies that there is some j € {1, . . . , m} such that xe{0) ^ for 



all € e §j . So, Lemma 7.4 insures that x{t) e K" for aU t e (0, T] 



Another special case of that of the feedback system (|l6|), and more generally when property ( [l9| ) 
holds. This property says that there is some jo S {1, . . . , to} such that, for every £ € §jg , it holds for each 
z e K"g that either (a) ze > 0, or (b) zi = an d gi{z) > (in which case also fl{z) — Pk{z)+gk{z) > 0). 



Applying with z = x{0) and using Lemma 7_A with this jo, we have: 



Corollary 7.6 Suppose that property ( fi9| ) holds, and consider any solution x : [0, T] R"q of 
Then, x{t) e «1 for aU < e (0, T]. " □ 

We also note, for further reference: 



Corollary 7.7 Consider a system (21), and pick any ^ G IR>o- Then, cither £ Eq or ^ belongs to 
some positive class. 



Proof. If ^ Eq, then Corollary 7.5, applied to a solution starting from ^, and forward invariance of 



classes, give that ^ is in the positive class containing x{T). 

8 Stability 

We start by establishing some useful estimates. 
Lemma 8.1 Define the following quadratic function: 



Q{Vi,---,Vm) ■■= ^^aij{r], - rijf . (51) 

1=1 3 = 1 



Then, there exists a constant k > such that 



Q((?i,...,g„) > n^^{qt- qjf (52) 

i=i j=i 



for aU (gi, . . . ,g„) G 



Proof. We first observe that 

Q{qi, ...,qm)=0 => qi = qm,i = i,---,rn-l. 

Indeed, obviously Q{qi, . . . , qm) = implies qi = qj for each pair i,j for which ^ 0. Now let J be the 
set of indices i such that qi = qm, and J its complement; as m G /, / 7^ 0. We need to see that J = 0. 
Suppose that J ^ %. The connectedness of the incidence graph of A provides an i G / and j £ J such 
that Oij ^ 0. Thus, qj = qi ^ qm, contradicting j G J. 

Let us introduce next a quadratic form in to — 1 variables: 

m—lm—l m—1 m—1 

i—1 j=l i—1 1 

which we denote as P{^i, . . . , Cm~i)- Since {rji — rjm) — {Vj ^ Vm) = Vi ^ Vj for all one has 

Qim^ ■ ■ • ,?7m) = ^'('71 - 77m, ■ • ■ jT/m-l ^ '7m) • 
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Note that P is positive definite: if P{qi, . . . , Qm-i) — 0, then Q{qi, . . . , Qm-i, 0) = 0, which as aheady 
observed implies that all qi = 0. Thus, there is some constant kq > such that 

m — 1 

P{pi, . . . ,Prn-l) > Kq 



pm— 1 



for all (pi, . . . jPm-i) ^ 1^"^ , which means that 

m— 1 

Q{qi, ■ ■ • , qm) > Ko ^ {qi - qrnf (53) 
1=1 

for all {qi, . . . , q,n) S R™. As {qi — qj)"^ < 2{qi — q^)^ + ^{qj — g,„)^ for all we may re-express the 
estimate ( ^3|) in the form (|5^), using a smaller constant k which depends only on kq and m. I 

The following estimate will be the basis of a Lyapunov function property to be established later. 

Lemma 8.2 There exist two continuous functions 

u : R" R" , c : R!;: ^ R+ 
such that, for every pair of points x, z in R" : 

{p{x) - ]{x)) < -c{z)S{x, z) + {v{p{x) - p{z)), f{z)) . (54) 

Proof. As B has full column rank, there is an to x n matrix B'^ (for instance, its pseudo-inverse) such 
that B*B = I. We let 

v{a) -.^ ((e<''^'->,...,e<^-->)i3#)' 

and 

Co(C) min e^''^''^^^^^ . 

j = l,...,m 

Now take any pair of positive vectors x, z. Denote, for each j = 1, . . . , m: 

qj := {bj , p{x) - p{z)) 

and observe that 

{bi,v{p{x) - p{z))) = e«* , « = TO 



so, using formula (24), 

m m 

{v{p{x) - p{z)), f{z)) =Y.T. "^.e^"^"'"^^^^^ i^" - e'O ■ (55) 
Therefore (writing g{x,z) — {v{p{x) — p{z)),f{z)) for simplicity): 

m m 

(p»~p(z),/(x)) = ^^a,,e<''-'^(-»(g,-9,) 

7n m 
7n m 

= E E "■^.e^"^"'''^'^^ (e'^ {q^ ~q,)~ (e"' - e«0) + (^, ^) (56) 

j=l 3 = 1 

in m 



^ -^EE«.,e<^-^^(^»(?.-'Z,)' +.9(^,^) (57) 
< -~2^a{z)^^a^j{q^-q.jf + g{x,z) 



1=1 j=i 

-^Co{z) Q{qi, ...,qm) + g{x, z) 
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where Q is the quadratic form in Lemma B.l. Equahty ( JSq ) foUows by adding and subtracting g{x, z) 
and using ([55|). To justify (p7[), we note first that, for each a > 0, the function M>o ^ M: 

fair) e'^(r _ a) - + e« + i(r - a)^ 

is always < (because /a(0) = -e"a — 1 + e° + (l/2)a^ < 0, fair) ~oo as r ^ +oo, and /^(r) = 
— e*" + (r — a) 7^ for aU r > 0). Now we use the inequality e"(r — a) — e*" + e° < —5(7" — a)^ in each 
term of the sum with a = qj and r = qi (recall 0^ 6^''^ '^*^^^^ > 0). 

Lemma B.l gives that Qiqi, . . . , qm) > K(5(a;, z). Thus, we may take c(z) := kco(2;)/2. I 



8.1 An Entropy Distance 

We will show the stability conclusions even if one does not impose the assumption that each 9i is onto 
R>o- Note that none of the results in Sections |^ and ^ used this condition. Thus, we suppose that each 
6'i : R [0, ai) (where < (7^ < 00) is locally Lipschitz, has 0^(0) = 0, satisfies \Piiy)\ dy < 00 (where 
Piiv) = lii0i(2/)), and its restriction to M>o is strictly increasing and onto [0, ai), with ai < +00. We let 
f5i := IntTi. (For the stability results, we will not even need to ask for condition (|28|).) 



For any fixed constant c G (—00, pi), and each i = 1, . . . , n, we consider the following function: 

Rii'r) / Piis)ds - cr . 



This function is a well-defined continuous mapping R>o M, continuously differentiable for r G (0,pi). 
Moreover, Rf achieves a global minimum at the unique r"^ = G R+ where Piir^) = c, decreases for 
r G [0, r'^], and increases to +00 for r > r'^. 

The following function will play a central role: 



W : R^o X M!^ ^ M : (a;, z) ^ RT^'^\^' 

i=l 



The above-mentioned properties of the functions R^^'"^''^ hiiply that 

x^z Wix,z) >Wiz,z) , (58) 

i.e., for each fixed z G M" , the function VF(-, z) has a unique global minimum, at z. Note also that the 
gradient of M^(-, z): 

dW 

— ix,z) = ipix)-piz))' (59) 

(defined for x G R") vanishes only at a; = z and that (since R^''''^'\xi) — > +00 if Xi +00) 

|a;| ^ +00 =^ Wix,z)^+oo (60) 
for every given z. As Wi-, z) is continuous, this implies that 

{x\Wix,z)<w} (61) 

is compact for every z and every w G R. 

Remark 8.3 In the standard setup pi = In, Wix, z) — J27=i 

Inz^. Then this formula, 

when states x are interpreted probabilistically in applications such as chemical networks, is suggested 
by "relative entropy" considerations. □ 
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8.2 Main Stability Results 

The next lemma will be applied later with S equal to the whole space or to a class, depending on the 
type of system. 

Lemma 8.4 Let S C M"q be a closed set, and pick x E S. Suppose that 

{p{x)-p{x),r{x)) < (62) 

is valid for all x £ 5p|M", x ^ x. Consider any ^ £ R"g for which the maximal solution x{i) of ( p^ ) 
with x(0) = ^ is included in S. 



1. If the system has the form ( |2l| ) then xit) is defined for alH > 0, and 

x(t) — > -EoUi^^} as t ^ +00 . 

2. If (^ holds, then x{t) is defined for aU i > 0, and 

x(t) ^ X as t +00 . 



Furthermore, in either case x is an equilibrium of (17), asymptotically stable relative to S. 

Proof. Fix 5*, x, and ^ as in the statement, and let x(-) be the maximal solution of ( p^ ) with x{0) = ^, 
defined a priori on some maximal interval [0,t*). We will use 

V{x) := W{x,x) -W{x,x) (63) 

as a Lyapunov-like function. By (p8|), this function is positive definite relative to the equilibrium x, 
i.e., V{x) > for all x £ K>05 ^^'^ ^(^) = if and only ii x = x. Moreover, V is proper, meaning 
that the sublevel sets {x \ V{x) < w} are compact, for all w £ M>o, by (|6l|). Finally, V is continuously 
differentiable in the interior M", and, using (p9|), 

VV{x)f*{x) < (64) 

whenever x £ 5 f] R+, 2: ^ 5, by (H). 

If the system has the form ( pl| ) and £, £ Eq, then ^ £ Eq (cf. Proposition 6.3), and thus, in that 
special case, x{t) is of course defined for all t > and converges to an equilibrium in Eq (it is constant, 
in fact). 

So, from now on, we sup pos e that either th e sy stem has the form ( ^T|) and ^ ^ Eo, or that ( |l9| ) holds. 
By, respectively, Corollary ^ and Corollary 7J, we know that x{t) £ S'p|R" for aU t £ [0,t*). So 
V{x{t)) is differentiable for t £ (0,t*), and dV{x{t))/dt < (by (64) if a;(i) ^ x, and obvious otherwise, 
since \7V{x) — 0), which means that V{x{t)) is nonincreasing. Since V is proper, this means that the 
maximal trajectory is precompact, and hence it is defined on the entire interval [0, +cxj), as claimed. 
Furthermore, the LaSalle Invariance Theorem implies that 

a;(t) r2+(^) as t +00 , 

where 51'*' (^) is the w-limit set of ^, which is a compact subset of {p \ V{p) = a}, for some a > 0. As the 
set S is closed, ^^^(0 ^ S. We pick any ( £ 51''' (^), and show that necessarily ^ = x or, in the case of 
systems (|l|), C ^ Eo- 
li C, — x, we are done, so we may assume from now on that ( ^ x. Similarly, if the system is of 
type (H^), we will assume that ( ^ Eq. We now derive a contradiction. 

Consider the forward trajectory z(t) starting from ^. Since r2''"(^) is invariant and a subset of S, 
z{t) £ S for all t, and z{t) ^ x for all t (since otherwise z{t) = x). Moreover, z{t) £ M" for all i > 
(using either Corollary |7.5| or Corollary [f^ ). Thus ( |64| ) says that dV{z{t))/dt < for all t > 0, which 
means that V{z{-)) is strictly decreasing. But this is a contradiction, since V{z{t)) = a. 

The stability statement is a simple consequence of the fact that F is a Lyapunov function (see e.g. 
section 5.7) relative to x for the dynamics restricted to SPlK". I 
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8.2.1 Proof of Theorem H 



Suppose r = n — m + 1, 7i, . . . , 7r are positive, ki, . . . ,kr G {1, . . . , n} are such that (14) holds, and j is 
such that (|l5|) holds. Fix any equilibrium x G 

We must prove that, for any given ^ G M"q, the maximal solution x{t) of ( [l6| ) with x{0) — is defined 
for all t > (the fact that the solution remains in R"q has already been established) and that x is a 
globally asymptotically stable equilibrium of (|l^). 

We apply Lemma B.4, with S = K'>o- Notice that 

(pix) ~ p{x),f{x)) < -c{x)5{x,x) < (65) 
for all X G M!f. by ( p^ (since f{x) — 0). The inequality is strict unless d{x,x) = 0, i.e., 

p{x) - p{x) G . (66) 

On the other hand, 



{p{x) - p{x),g{x)) = ^'^i{xk, - Xke) {p{x) - p{x),ek,) 

r 

= '^leixki - Xkf) (pkeixke) - Pkeixke)) < 

e=i 

where the last inequality follows from the fact that each pi is an increasing function, and the inequality 
is strict unless 

{p{x) - p{x),eki) Pktixke) - Pkeixkt) = (67) 
for £ = I . . .r. Thus {p{x) — p{x), f*{x)) equals 

{pix) ~ pix)Jix)) + {pix)~ p{x),g{x)) < 0, 



and this inner product can only vanish if both (|66|) and (|6^) hold, which implies, because of (^^, that 
p{x) — p{x), i. e., X — X. This means that ( |6^ holds. As the system (|l^ ) is a system ( |l7|) for which ( [l9| ) 
holds. Lemma 8.4 provides the proof of global stability. I 



8.2.2 Proof of Theorem g 



Eq — Eq, 



Consider any system (21) and any maximal solution x{-). By Corollary 7.7, either x{t) = G 
in which case of course x{t) — *■ E, or there is a positive class S such that x{t) G S for all t, which we 
assume from now on. Let x = xg- Notice that ( |65| ) again holds, for all x G M" . Let us now specialize to 
systems (^. We claim that 

{p{x) ^ p{x)J{x)) < (69) 

(i.e., property ( |6^ ) for (||)) is vahd for all x G S'P|R" , x ^ x. Since c{x) > 0, by ( |65| ) all we need to show 
is that S{x, x) ^ 0. But the only way that (5(x, x) could vanish is if x is an equilibrium (cf. Corollary 5.3), 
and uniqueness of equilibria in 5p|M!f: then gives x = x. So, Part 1 in Lemma B.4 gives x(t) —> E. I 



8.2.3 Proof of Theorem g 

Let be a positive class, and pick x — xs- Define: 

6s{x) := ^cix)"^ 5{x,x) 

restricted to x G iSpR". Since c{x) > and S{x,x) = if and only if x is an equilibrium, which by 
uniqueness means x ~ a;, we have Ss{x) > for all x ^ x. take now any collection of functions {Ai^} 
such that (|2|) holds. 
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Notice that (|65|) again holds, for all x G R". On the other hand, for x G S'p|] 



{p{x) - p{x),g{x)) = Ay(x) (p(a;) -p(a;),6, 

i=l 3 = 1 

< ^ySs{x)^/^{x,x) 

< -c{x)S{x,x), 

where the first inequality follows from the Cauchy-Schwartz inequality, the estimate (|2^), and the defi- 
nition of 6. Therefore {p{x) — p{x), f*{x)) equals: 

{pix) - p{x), f{x)) + {p{x) - p{x), g{x)) < -i c{x)5{x, x) (70) 

for all x G S'PlM", and this expression is negative when x ^ x. Thus, hypothesis ( |6^ ) in Lemma 8.4 
holds, with S as given. Part 3 (asymptotic stability of x) follows from the Lemma. 

We pick any G S*, and consider the ensuing maximal solution. By Part 1 of Lemma the solution 
is defined for all t > (proving Part 2 in Theorem ^ and x{t) Eq[J{x} as < ^ +oo. Since Eg and 
{x} are at positive distance, this means that either x{t) — > Eq or x{t) — > x. In the first case, 5* being 
closed implies that x{t) Eof]S. Thus, if Eq — Eq does not intersect 5', the only possibility is that 
x{t) X. Conversely, if S'P] -Eq ^ then xs is not globally asymptotically stable relative to S. This is 
clear, since ii ( € Eof] S ^ (l> then ^ (being an equilibrium) is not attracted to x. This proves Part 4 of 
Theorem |. I 

Observe d{x,x) = S{x,Tr{x)) is of class C'^ if hypothesis (H^) holds. One may pick easily a smooth 
lower bound for the function c, and in that manner obtain a class C*^ function Ss- 



8.2.4 Proof of Parts b and c in Theorem |^ and Lemma 2.2 



Pick any positive class S. We let Ss be any function as in the statement of Theorem ^; as (x) = 0, 
the hypotheses of that theorem are verified. In this manner, all conclusions in Theorem ^ as well as 



Lemma 2.2 are established. 



Remark 8.5 A different approach to the proof of Theorem not using LaSalle invariance, can be 
based upon the fact that V{x{t)) can be shown to decrease strictly along every trajectory (not merely 
in the positive orthant), and this proof applies as well in the multiple- linkage class case described in 
Remark ^.5[ (The fact that all trajectories approach equilibria, in the multiple-linkage class case, may 
also be proved directly. The LaSalle argument proceeds in basically the same manner; the critical step is 
to study those trajectories that remain in a level set of the Lyapunov function and start on the boundary. 
In the case of a single class, such trajectories are either equilibria or enter immediately the interior M" 
(Corollary |7.5| ). In the general case, there may be a block which is "turned off", that is, fi{x{t)) = 
for some i G {1, . . . , L}. In that event, we are reduced to studying a system with L — 1 classes, and the 
same Lyapunov function (restricted to a possibly smaller state space, after dropping those state variables 
which only appear in the "off" reaction) can be used inductively. □ 



8.3 A Remark on Exponential Stability 

Suppose that hypothesis (Hi) holds. Recalling the definition ( |6^ ) of the Lyapunov function V, and 
the form of its gradient given by (^9|), {p{x) — p{x))' , we know that the Hessian of ^ at a; is given by 
Q = diag(pi(a;i), . . . , p'^{xn)) and is therefore nonsingular. Since the gradient of V vanishes at a; = a;, we 
have a Taylor expansion V{x) ~ [x~ x)'Q{x ~ x) + o{\x ~ x\ ), and thus, for all x in some neighborhood 
of X, ci \x — x'f' < V{x) < C2\x — x\^ for some positive constants ci and C2. If we let d2 — d2(x)/2 
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in (|3^), we also have that S{x, x) > d2\x — x\^ for aU x in a neighborhood of x. So inequahty ( [70| ) gives 
us that, along trajectories in the class containing x, as long as x{t) is near x we have 

^ffi^ = {p{x)-p{xirix)) < -cVixit)), 

where c = c{x)d2Ci/2. Integrating, V{x{t)) < e~''*V{x{0)), so we obtain an estimate \x{t) — x\'^ < 
(c2/ci)e~'^* |a;(0) — x\^ for all trajectories of system which start near x. In other words, (relative) 
stability is in fact exponential. I 
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